Pattern selection in a lattice of pulse-coupled oscillators 



X. Guardiolala and A. Diaz-Guileral] 
Departament de Fisica Fonamental, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain 

(February 1, 2008) 









a 

o 
o 



> 

m 

o 



I 

O 

o 



X 



We study spatio-temporal pattern formation in a ring of 
A'' oscillators with inhibitory unidirectional pulselike interac- 
tions. The attractors of the dynamics are limit cycles where 
each oscillator fires once and only once. Since some of these 
limit cycles lead tothe same pattern, we introduce the con- 
cept of pattern degeneracy to take it into account. Moreover, 
we give a qualitative estimation of the volume of the basin 
of attraction of each pattern by means of some probabilistic 
arguments and pattern degeneracy, and show how are they 
modified as we change the value of the coupling strength. In 
the limit of small coupling, our estimative formula gives a 
perfect agreement with numerical simulations. 
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I. INTRODUCTION 

The study of the collective behavior of populations of 
interacting nonlinear oscillators has attracted the inter- 
est of physicists and mathematicians for many years since 
they can be used to modelize several chemical, biologi- 
cal and physical systems 00]. Among them, we should 
mention cardiac pacemakers cells [||, integrate and fire 
neurons B and other systems made of excitable units 
g. Most of the theoretical papers that have appeared 
in the scientific literature deal with oscillators interact- 
ing through continuous-time couplings, allowing them 
to describe the system by means of coupled differential 
equations and apply most of the modern nonlinear dy- 
namics techniques. More challenging from a theoretical 
point of view is to consider a pulse-coupling, or in other 
words, oscillators coupled through instantaneous inter- 
actions that take place in a very specific moment of its 
period. The richness of behavior of these pulse-coupled 
oscillatory systems includes synchronization phenomena 
pI , spatio-temporal pattern formation M] (we could men- 
tion, for instance, traveling waves H, chessboard struc- 
tures [0, and periodic waves [nO[ ), rhythm anihilation 
||11| , self-organized criticality [||!7. . 

Most of the work on pattern formation has been done 
in mean-field models or populations of just a few oscilla- 
tors. However, such restrictions do not allow to consider 
the effect of certain variables whose effect can be cru- 
cial for realistic systems. The specific topology of the 
connections or geometry of the system are some typi- 
cal examples which usually induce important changes in 
the collective behavior of these models. Pattern forma- 
tion usually takes place when oscillatory units interact in 
an inhibitory way, although it has also been shown that 



the shape of the interacting pulse, when the spike lasts 
for a certain amount of time, or time delays in the in- 
teractions can lead to spatio-temporal pattern formation 
also in the case of excitatory couplings |lj,|l^. Only 
recently, general solutions for the general case, where 
the patterns existence and stability is proved, have been 
worked out [|l2|,^. The aim of this paper is to study 
some pattern properties and get a quantitative estima- 
tion of the probability of pattern selection under arbi- 
trary initial conditions or, in the language of dynamical 
systems, the volume of the basin of attraction of each 
pattern. Keeping this goal in mind, we will use the gen- 
eral results given in |13(| where assuming a system de- 
fined on a ring the authors developed a mathematical 
formalism powerful enough to get analytic information 
of the system. Not only about the mechanisms which are 
responsible for synchronization and formation of spatio- 
temporal structures, but also, as a complement, to proof 
under which conditions they are stable solutions of the 
dynamical equations. 

Despite the apparent simplicity of the model, some ring 
lattices of pulse-coupled oscillators are currently used to 
modelize certain types of cardiac arhythmias where there 
is an abnormally rapid heartbeat whose period is set by 
the time that an excitation takes to travel the circuit pGl . 
Moreover, there are experiments where rings of a few R15 
neurons irora Aplysia are constructed and stable patterns 
are reported [0 . Our Id model allows us to study ana- 
lytically the most simple patterns and understand their 
mechanisms of selection. 

The structure of this paper is as follows. In Sec II we 
review the model introduced in [Of as well as set the 
notation used throughout the paper. In Section III we 
study some pattern properties which will be useful for, 
in Section IV, propose an estimation of the probability of 
selection of each pattern. In the last section we present 
our conclusions. 



II. THE MODEL 

Our system consists in a ring of (A^-|- 1) pulse-coupled 
oscillators. The phase of each oscillator (j)i evolves lin- 
early in time 



dt 



1 



Vi = 0, 



,iV 



(1) 



until one of them reaches the threshold value (jjth = 1- 
When this happens the oscillator fires and changes the 
state of its rightmost neighbor according to 



> 1 => 



Ui+l 



fl(f), 



'i+1 



subjected to periodic boundary conditions, i.e. A^+1 = 0, 
and where e denotes the strength of the couphng and 
fi — 1 + £. Where we have assumed that, from an ef- 
fective point of view, the pulse-interaction between os- 
cihators, as well as the state of each unit of the system, 
can be described in terms of changes in the phase, or 
in other words, in terms of the so called phase response 
curve (PRC), £0 in our case. A PRC for a given oscil- 
lator represents the phase advance or delay as a result 
of receiving an external stimuli (the pulse) at different 
moments in the cycle of the oscillator. We will assume 
e < througout the paper, as we are only interested 
in spatio-temporal pattern formation and e > always 
leads to the globally synchronized state |jl3]. This linear 
PRC has physical sense in some situations. For instance, 
it shows up when we expand the non-linear PRC for the 
Peskin model of pacemaker cardiac cells Q in powers of 
the convexity of the driving or in neuronal modelling mm . 
In practice, however, this condition can be relaxed since 
a nonlinear PRC does not change the qualitative behav- 
ior of the model provided the number of fixed points of 
the dynamics is not altered. Moreover, a linear PRC has 
the advantage of making the system tractable from an 
analytical point of view. 

Let us describe the notation used in the paper. The 
population is ordered according to the following criterion: 
The oscillator which fires will be always labeled as unit 
and the rest of the population will be ordered from 
this unit clockwise. After the firing, the system is driven 
until another oscillator reaches the threshold. Then, we 
relabel the units such that the oscillator at = 1 is now 
unit number 0, and so on. This firing -I- driving (FD) 
process for A^ -f- 1 oscillators can be described through a 
suitable transformation 



0' = Tfc(0) = r+Mfe0, 



(2) 



where M^ is a A^ x A^ matrix, is a vector with A^ compo- 
nents, 1 is a vector with all its components equal to one 
and k stands for the index of the oscillator which will fire 
next. We call this kind of transformation a firing map, 
and we have to define as many firing maps as oscillators 
could fire, that is, index k must run from k = 1 (0i fires) 
to A^ (0JV fires). For example, in the N + 1 — 4 oscilla- 
tors case we have that the firing map corresponding to 
the FD process where 02 is the next oscillator which do 
fire, 
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and so on. Once we have defined all possible firing maps 
for a given number of oscillators we can proceed to deal 
with the attractors or fixed points of the system dynam- 
ics. As has been proved in |n3[ these fixed points must be 
cycles of A^ -I- 1 firings. We define a cycle as a sequence 
of consecutive firings where each oscillator fires once and 
only once. Mathematically, each cycle is described by 
means of a return map. The return map is the transfor- 
mation that gives the evolution of during a cycle and 
is the composition of all firing maps involved in the firing 
sequence of that cycle 



b' = Te, o T,, 



or,„^,(0)= R, + ] 



(4) 



where Tc- o Tc (0) is the usual composition operation 
T,,(Te^.(0)) and 



Rr = l 



Cn 



iU 



i=ci j=ci 



and 



Cn + 1 



Note that not all possible combinations of firing maps are 
allowed, just those ones whose indices q sum p(N + I) 
without any partial sum equal to q{N + 1), where p > q 
are positive integers. 

As all firing maps are linear transformations, return 
maps are also linear. There are A! possible cycles in the 
A+1 oscillators case (all permutations of firing sequences 
with the initial firing oscillator 0o fixed). Following our 
previous example, for the four oscillators case all possible 
firing sequences and their associated return maps are 



A: 0,1,2,3 
B : 0,1,3,2 
C: 0,2, 1,3 
£>:0,2,3,1 
£;:0,3,1,2 
F: 0,3, 2,1 



Ti o Ti o Ti o Ti 
T20T30T20 Ti 
T10T20T30 T2 
T30T20T10 T2 
T20T10T20 T3 
T30T30T30 T3 



Now, in order to find the attractors of the dynamics, we 
must solve the fixed point equation 



0: = i?e+Me0:, 

for every cycle c. Formally, 

$:=Rc-{I-M,)-\ 



(5) 



(6) 



As was shown in jl^, there are A^ different stable solu- 
tions to the whole set of fixed point equations. Their 
stability is assured by the fact that e < 0, since it guar- 
antees that all eigenvalues of Mc lie inside the unit circle 
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FIG. 1. Graphic representation of the patterns solution 
of Eq. (7) for a small value of the coupling strength |e| at 
the beginning of the cycle (we must keep in mind that spa- 
tio-temporal patterns are dynamical structures that evolve in 
time). The leftmost square represents 0o and the rightmost 
4>N, and their phase is visualized in a greyscale where black 
means th — 1 and white d> = 0. 



for all cycles c. In our four oscillators example these so- 
lutions are 



0A-(1,- 



:) 



:) 



4+3£' 4+3£' 4+3e^ 

<f*B ^ <Pc ^ <P*D = <P*E = (1' 2Ti'^'2Ti) 

rh* - n ^- 2+£ 3+£] 
^F V-^i 4+ei 4+e: 4,+e > 

Which are a kind of four-oscillators traveling wave, chess- 
board and inverse traveling wave structures. 

From now on we will label such solutions with index 
TO (to = 1...N) since their first component always satisfy 

^^ N + l + me' ^ ' 

Therefore, in the example, we relabel patterns (j)\ as ttt- = 



"Bi v^c vd^ ve 



as m — 2 and 



as TO, = 1. 



Since there are A^! possible cycles and N solutions to 
Eq. (7) there will be some fixed points or patterns which 
will appear more than once, so, we shall use C{N -I- 1, to) 
to characterize these degeneracies. In the example, the 
values of the degeneracies are C(4, 1) = C(4,3) = 1 and 
C(4, 2) = 4. In general, patterns which are solutions of 
cycle consisting in the iterative application of the same 
firing map (like A and F in our example) have no period- 
icities whereas the ones solution of mixtures of differents 
firing maps (B,C,D and E) have some periodic structure 
that are also solution of Eq. (7) for a case with less os- 
cillators. In Fig. 1 we can visualize the solutions for 
A^ + 1 = 2,3 and 4 oscillators and realize that solution 
TO = 2 for the four oscillators case is a periodic composi- 
tion of solution m — 1 for the two oscillators case. 



III. PATTERN PROPERTIES 

As we have seen, the stability of all patterns solution 
of Eq. (6) is guaranteed by the fact that e < 0, but 



the existence of such solutions is not ensured. In fact, 
for small values of the coupling strength |e| all patterns 
do exist, but, as we increase it, some patterns disappear. 
The reason is that the solution loses its physical meaning 
because 0| > 1. Their first component is always the one 
that becomes larger than unity earlier and this happens, 
for each to, and according to Eq. (9), when 



e < £„ 



1 



A^ + 1 



(8) 



Our coupling strength range of interest ends at e = — 1 , 
since at e < — 1 we always find the same pathological 
dynamics which does not have any physical or biological 
sense. Realistic couplings never reach such higher values. 
Therefore, as e runs from to —1, all patterns whose m 
satisfy m > ^^^^^, disappear. 

There is another interesting pattern property which 
has to do with the calculation of the pattern degeneracy 
C{N -I- IjTO.). In principle, to calculate such degenera- 
tion, we should solve fixed point Eq. (6) for all possible 
cycles and count how many of them lead to the same pat- 
tern. Although for few oscillators the problem is quite 
straightforward, as we deal with higher and higher num- 
ber of oscillators, the number of cycles increases (it grows 
as Nl) and solving Eq. (6) becomes more difficult. For- 
tunately, there is another way of calculating C(A^ + 1, to,) 
which reduces the problem to a combinatorial question. 
Lets show it through an example, in the previous four 
oscillators case, if we count, for each firing sequence, the 
number of oscillators which have received the pulse be- 
fore firing, we can easily realize that this number is the 
same as its value of m 
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0,1,2,3 


TO 


= 3 
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0,1,3,2 


TO, 


= 2 
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0,2,T,3 


TO 


= 2 
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0,2,3,T 


TO 


= 2 
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0,3,T,2 


TO 


= 2 


F 


0,3,2, T 


TO 


= 1 



Here an upper bar means that the oscillator has already 
received a pulse during the cycle. The point is that 
it turns out that every pattern to, corresponds to a se- 
quences of firings involving exactly to oscillators that, 
when they do fire, had already received a pulse from their 
leftmost neighbor. Therefore, this property (we have 
checked for several values of TV -f 1) allows us to associate 
every cycle with the pattern it leads to, just by count- 
ing these kind of firings. Now, calculating C{N + 1,171) 
becomes a straightforward matter. In Table I we have 
computed C{N -I- 1, to) for several values of iV -I- 1. 

Apart from brute force counting, degeneracy distribu- 
tion C{N + l,m) can also be determined from the fol- 
lowing relation 

C{N + 1, to) == mC{N, to) + 

(iV + l-m)C(A^,m-l), (9) 



TABLE I. Pattern degeneracy C{N + l,m). First column 
stands for the number A'' + 1 of oscillators and first row for 
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14608 


88234 


156190 


88234 


14608 


502 


1 



for 2 < m < N — \. This recursion relation is closed by 
C{N + 1,1) ^C{N + 1,N)^ I, (10) 

which correspond to the firing sequences 



0,7V, (iV- 1)... 2,1 and 0, 1, 2...(iV - l),iV, 

respectively. 

From the previous relations one can deduce by induc- 
tion the symmetry of the distribution with respect to its 
extremes at tti = 1 and m = N 



and 



C(A^ + 1,to) = C(iV + l,iV+l-TO), (11) 



^C(A^+l,m) = A^!. (12) 



Another interesting property is the period A^"*"^ of each 
spatio-temporal pattern m. Since all oscillators are in 
a phase-locked state, they must oscillate with the same 
period. Then, as the intrinsic period of each oscillator is 
one, and when any oscillator receives the delaying pulse 
from its neighbor it has a phase equal to 0^, one can 
easily realize that the effective period is 



A 



N+l 



1 



A^ + 1 + 2me 
A^ + 1 + TO£ ■ 



(13) 



Therefore, the larger the value of m, the longer the period 
of its associated pattern. It is important to notice that we 
have not fixed the value of such periods (each pattern has 
its own which is different from the others), since there are 
some authors who fix all periods equal to some constant, 
and use it as a condition to find the structures [17|. 



IV. PATTERN SELECTION 

Once we have characterized all spatio-temporal pat- 
terns, we proceed to find some general formula which give 
us some estimation of the probability of each pattern to 



be selected, or in other words, an estimation of the vol- 
ume of its basin of attraction. In order to achieve this ob- 
jective, we should understand the mechanism which lead 
to the selection of a certain spatio-temporal structure and 
how is it modified as the parameters of the model (e in 
our case) change. 

There is an easy and straightforward way to get the 
essential features of this mechanism assuming that the 
probability of one oscillator to fire next is, basically, pro- 
portional to its phase (that is, if it has a phase slightly 
below 1 it has a higher probability to be the next fir- 
ing oscillator, whereas if it has a smaller phase, it will 
rarely fire next). Imagine the phases of all oscillators 
randomly distributed over the interval (0,1). Then we 
let the system evolve till one of the oscillators reaches a 
phase (pi = 1 and emits a pulse that is received by its 
rightmost neighbor which lows its phase by an amount 
e0i+i. Now we assume that all phases are again ran- 
domly distributed over (0,1) except the one which re- 
ceived the pulse whose phase is distributed over (0, 1-f-e). 
So, we get rid of memory effects (we know the oscillator 
that has fired should, now, have a phase equal to zero) 
and just keep in mind if each oscillator has received a 
pulse or has not. Therefore, the point is that under this 
conditions, the probability that one oscillator which has 
still not received a pulse do fire is some constant and, 
on the other hand, for the ones which had, is this con- 
stant times the factor (1 -I- e). Then, we can character- 
ize the probability of having some cycle just by recalling 
how many oscillators do fire having previously received 
a pulse during that cycle. Basically, this probability is 
proportional to (1-l-e)" where n stands for the number of 
oscillators which do fire having already received a pulse 
(the product of all constant terms will be absorbed in a 
normalization factor). This approach, where we assume 
all firings as almost-independent events, can be viewed 
as a kind of mean-field approximation. Then, as has been 
shown before, since cycles leading to the same pattern 
m always exactly have m oscillators that do fire having 
received the interacting pulse, we can give an estimation 
of the probability for pattern m selection in the A^ -I- 1 
oscillators case 



pf^+\e)c^N{e)C{N + \,m){\ + eY 



(14) 



Here M{e) is chosen so that summation of the probabili- 
ties over m gives 1 



l^Pm 



N+l(^\ - 



(e) = 1. 



(15) 



In the limit of small coupling strength £ — > 0, which 
is the more interesting case for the majority of physical 
and biological systems, one can assume that interaction 
plays almost no role when pattern selection takes place. 
That is, the fact that one oscillator has received the pulse 
from its neighbor does not low its probability to fire as 
the pulse does not modify appreciably its phase. Then, 
we can consider that all cycles have approximately the 



same probability to be selected, (1 + e)™ -^ 1, and only 
pattern degeneracy has to be considered to get a good 
estimation of p^"*"^ 



P\ 



m 



(16) 



The dominant pattern, that is, the one which has the 
larger probability to be selected coincides with the mean 
value of m (due to the symmetric behavior of C{N + 
l,m)). 



< TO >jv+l=^l 



C(A^ + 1,to) _ A^+1 



(17) 



For an odd number of oscillators < to, >n+i docs not 
exist and we have a competition between the two closest 
patterns to = N/2 and to = (TV + 2)/2. Recall that 
the most probable patterns turn out to be the ones with 
" shortest wavelengths" , a fact that was already reported 
in simulations of these sort of systems Q. In Figs. 2 and 
3 we check this new approximation for the A^ + 1 = 10 
and 9 case and realize that expected results are in good 
agreement with simulations data. 

There also is the interesting question of how does this 
probability distribution modifies when the number of os- 
cillators increases. In Fig. 4 we show p,„ for different 
values of TV + 1. Since there are more possible values of 
TO. available, as we increase N +1, p^"*"^ diminishes. The 
distribution also gets narrower as we increase A^ + 1 and 
this becomes clear when one studies the variance of Pm- 
It can be found that 



< m^ >Ar+i = y^ i 



,C(JV+l,m) 

iV! 



(N + l)^ N+l 
4 ^ 12 



(18) 



We could not prove this without an explicit expression 
for C{N + 1, m) but we have checked it N up to 170. 
Therefore 



2 



N+l < m >N+i 



12 



6 



(19) 



It turns out that for a large number of oscillators almost 
all initial conditions lead to a pattern whose m approx- 
imately falls in the interval < m >n+i ±^/<l7^>Jv+T■ 
In order to compare it for different number of oscillators 
we have to normalize m dividing by A^ + 1. In that case, 
one observes that (t%^-^ ^ ^/VN + 1 so that as we in- 
crease A^ -I- 1, the spread of p^+^ diminishes getting the 
distribution sharpened. 

As Eq. (14) does not take into account the disappear- 
ance of the different patterns to at the different values of 
e^ predicted by Eq. (9), it can not give a good quanti- 
tative estimation of pattern selection for higher coupling 
values. Nevertheless we can expand Eq. (14) to the lead- 
ing order in e. For small e, p^^^ ^'"^ approximated by 




m 

FIG. 2. Estimation oipm for small coupling strength in the 
case of an even number of oscillators A*' 4- 1 = 10 and for a 
values of e equal to —0.01 and —0.1. Dashed line follows the 
theoretical values estimated through Eq. (16). We can real- 
ize that the smaller |£| is, the more accurate our estimations 
are. The most probable pattern is m = (A'^ -I- l)/2 and the 
probability for the patterns near the extremes is almost zero 
due to the fast decay oi pm there. 
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FIG. 3. The same as Fig. 2 but now for an odd number 
of oscillators A'^ -I- 1 = 9. We can realize that there is not a 
peak anymore, instead, almost all probability is concentrated 
in the two competing patterns m — N/2 and m — {N + 2)/2. 
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N+l 



C{N + l,m) 



(1 + (to 



N + l 



e). 



(20) 



In Fig. 5 we compare this approximation with simulated 
data. The slopes near e = do agree with Eq. (21). 
In our simulations we calculate the probability of each 
pattern to be selected just by counting how many real- 
izations (with ipQ = 1 and the rest of oscillators with ran- 
dom initial conditions) lead to each pattern to. and divide 
over the total number of realizations. Although we only 
have a good quantitative estimation of p^^^ for small 
values of e, Eq. (15) catches the two basic mechanisms 
responsible of pattern selection. On the one hand, it is 
clear that for higher values of the coupling strength \e\, 
when one oscillator receives a pulse, it lows its phase to 
almost zero and, consequently, its firing probability also 
does. Therefore pattern selection probability p^+^(e) is 
strongly controlled by the number of oscillators which 
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FIG. 4. pm calculated by means of Eq. (17) for different 
values of A''+ 1. The distribution gets narrower and the height 
of its peak diminishes as we increase the number of oscillators. 
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FIG. 5. Estimation of p^^^ for a system of N+l=6 oscil- 
lators by means of the linearized Eq. (21). For small e the 
slopes match with simulated data. 



have to fire having already received a pulse, that is, the 
probabihstic factor (1 + e)™. As a consequence, p^'^^ 
begin to decrease sooner vi^hen \e\ increases, the larger m 
is. On the other hand, for small values of the coupling 
strength, interaction plays almost no role and p^'^^{£) is 
dominated by the degeneracy factor C{N + 1, m). There- 
fore p^^^(£) for the different values of m are basically 
ordered as C{N + \,m). In Fig. 6, 7 and 8 we show 
results from simulations of p^^'^(e) for different number 
of oscillators. 



V. CONCLUSIONS 

In this paper we have studied some properties of the 
spatio-temporal patterns that appear in a ring of pulse- 
coupled oscillators with inhibitory interactions. We have 
focused our attention in estimating the probability of se- 
lecting a certain pattern under arbitrary initial condi- 
tions and have shown the two basic mechanisms respon- 
sible of that: the degeneracy distribution C{N + l,m), 
for small values of e, and m, the number of oscillators 
that do fire having already received a pulse, for higher 
values of e. According to this, the different probabilities 
of selecting pattern m start being distributed following 
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FIG. 6. Pm(e) in the case of A^ -I- 1 = 4 oscillators. In 
this case, we see that although for small coupling chessboard 
(m = 2) is the dominant pattern, the inverse traveling wave 
(m = 1) is the most probable pattern for higher values of the 
coupling strength. Simulations are done over 2500 realiza- 
tions. 
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FIG. 7. The same as Fig. 6 for iV + 1 = 5. 
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the degeneracy distribution C{N + l,m), and, as e de- 
creases, these probabihties diminish in a hierarchical way: 
the larger the value of m, the sooner its selection prob- 
ability is going to decrease, so that only patterns with 
smaller m will survive for higher values of e. Moreover, 
some of the structures disappear, at the different values 
of £^, during this process. We have found out an approx- 
imation formula for p^^^^i^) which takes into account all 
these mechanisms and gives us a quantitative estimation 
of the different selection probabilities for small e. 

The estimation of the volume of the basin of attraction 
of each spatio-temporal pattern m also gives us an idea 
of the stability of the different structures with respect to 
additive noise fluctuations (for instance, we can add some 
random quantity -q to all phases after each firing event 
or a continuous-time rj{t) in the driving). Simulations 
of arrays of noisy pulse coupled oscillators showed that 
our most probable patterns were also the most stable . 
The present paper only concerns spatio-temporal pattern 
formation in a ring of oscillators, nevertheless, all results 
are trivially generalized to bidirectional couplings. Al- 
though the question of what happens when dealing with 
higher dimension lattices remains opened, some simula- 
tions results in 2d showed that almost all realizations 
lead to a chessboard pattern in analogy with our results 
in the ring. That makes us believe we have caught the 
basic features of the problem in our Id model. 
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